Asymptotic Approximations to Truncation 
Errors of Series Representations for Special 
Functions 

Ernst Joachim Weniger 

Institut fiir Physikalische und Theoretische Chemie 
Universitat Rcgcnsburg, D-93040 Regensburg, Germany 
j oachim . wenigerOchemie . uni-regensburg . de 

Summary. Asymptotic approximations (n — > oo) to the truncation errors r„ = 
— X^J^o '^'^ infinite series X^J^o '^^ special functions are constructed by solving 
a system of linear equations. The linear equations follow from an approximative 
solution of the inhomogeneous difference equation Zir„ = a„+i. In the case of the 
remainder of the Dirichlet series for the Riemann zeta function, the linear equations 
can be solved in closed form, reproducing the corresponding Euler-Maclaurin for- 
mula. In the case of the other series considered - the Gaussian hypergeometric series 
2-Fi(a, b; c; z) and the divergent asymptotic inverse power series for the exponential 
integral Ei (z) - the corresponding linear equations are solved symbolically with the 
help of Maple. The practical usefulness of the new formalism is demonstrated by 
some numerical examples. 

1 Introduction 

A large part of special function theory had been developed already in the 
19th century. Thus, it is tempting to believe that our knowledge about special 
functions is essentially complete and that no significant new developments are 
to be expected. However, up to the middle of the 20th century, research on 
special functions had emphasized analytical results, whereas the efficient and 
reliable evaluation of most special functions had been - and to some extend 
still is - a more or less unsolved problem. 

Due to the impact of computers on mathematics, the situation has changed 
substantially. We witness a revival of interest in special functions. The general 
availability of electronic computers in combination with the development of 
powerful computer algebra systems like Maple or Mathematica opened up 
many new applications, and it also created a great demand for efficient and 
reliable computational schemes (see for example or Section 13] 

and references therein). 



2 Ernst Joachim Weniger 

Most special functions are defined via infinite series. Examples are the 
Dirichlet series for the Riemann zeta function, 

oo 

Cis) = ^(^ + f)-^ (1) 

which converges for Re(s) > 1, or the Gaussian hypergeometric series 

2F,ia,b;c;z) = i^M^ z% (2) 

which converges for \z\ < 1 . 

The definition of special functions via infinite series is to some extent highly 
advantageous since it greatly facilitates analytical manipulations. However, 
from a purely numerical point of view, infinite series representations are at 
best a mixed blessing. For example, the Dirichlet series (Q) converges for 
Re(s) > 1, but is notorious for extremely slow convergence if Re(s) is only 
slightly larger than one. Similarly, the Gaussian hypergeometric series (O 
converges only for \z\ < 1, but the corresponding Gaussian hypergeometric 
function is a multivalued function defined in the whole complex plane with 
branch points at z = 1 and oo. A different computational problem occurs in 
the case of the asymptotic series for the exponential integral: 

oo 

ze-Si(z) ^ ^(-l/z)™m! = 2^0(1, 1; , z ^ 00 . (3) 

m— 

This series is probably the most simple example of a large class of series that 
diverge for every finite argument z and that are only asymptotic in the sense 
of Poincare as z — > 00. In contrast, the exponential integral i?i(z), which has 
a cut along the negative real axis, is defined in the whole complex plane. 

Problems with slow convergence or divergence were encountered already in 
the early days of calculus. Thus, numerical techniques for the acceleration of 
convergence or the summation of divergent series are almost as old as calculus. 
According to Knopp [SJ p. 249], the first systematic work in this direction can 
be found in Stirling's book IT7|, which was published already in 1730 (recently, 
Tweddle .20^ published a new annotated translation), and in 1755 Euler jHI 
published the series transformation which now bears his name. 

The convergence and divergence problems mentioned above can be for- 
malized as follows: Let us assume that the partial sums s„ = X]"=o '^'^ ^ 
convergent or divergent but summable series form a sequence {sn}^=o whose 
elements can be partitioned into a (generalized) limit s and a remainder or 
truncation error r„ according to 

Sn = s + rn, n e No . (4) 



This implies 
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r„ = - ^ , n G No . (5) 

v—n+l 

At least in principle, a convergent infinite series can be evaluated by adding 
up the terms successively until the remainders become negligible. This ap- 
proach has two obvious shortcomings. Firstly, convergence can be so slow that 
it is uneconomical or practically impossible to achieve sufficient accuracy. Sec- 
ondly, this approach does not work in the case of a divergent but summable 
series because increasing the index n normally only aggravates divergence. 

As a principal alternative, we can try to compute a sufficiently accurate 
approximation f„ to the truncation error r„. If this is possible, f„ can be 
eliminated from s„, yielding a (much) better approximation s„ — f„ to the 
(generalized) limit s than s„ itself. 

This approach looks very appealing since it is in principle remarkably 
powerful. In addition, it can avoid the troublesome asymptotic regime of large 
indices n, and it also works in the case of divergent but summable sequences 
and series. Unfortunately, it is by no means easy to obtain sufficiently accurate 
approximations f„ to truncation errors r„. The straightforward computation 
of r„ by adding up the terms does not gain anything. 

The Euler-Maclaurin formula, which is discussed in Section |21 is a prin- 
cipal analytical tool that produces asymptotic approximations to truncation 
errors of monotone series in terms of integrals plus correction terms. Unfortu- 
nately, it is not always possible to apply the Euler-Maclaurin formula. Given 
a reasonably well behaved integrand, it is straightforward to compute a sum 
of integrand values plus derivatives of the integrand. But for a given series 
term a„, it may be prohibitively difficult to differentiate and integrate it with 
respect to the index n. 

In Section 13 an alternative approach for the construction of asymptotic 
approximations {n —^ oo) to the truncation errors r„ of infinite series is pro- 
posed that is based on the solution of a system of linear equations. The linear 
equations exist under very mild conditions: It is only necessary that the ratio 
an+2/o.n+i or similar ratios of series terms possesses an asymptotic expansion 
in terms of inverse powers l/{n -j- a) with a > 0. Moreover, it is also fairly to 
solve these linear equations since they have a triangular structure. 

The asymptotic nature of the approximants makes it difficult to use them 
also for small indices n, although this would be highly desirable. In Section^ it 
is mentioned briefly that factorial series and Fade approximants can be helpful 
in this respect since they can accomplish a numerical analytic continuation. 

In Section [51 the formalism proposed in this article is applied to the trun- 
cation error of the Dirichlct series for the Riemann zeta function. It is shown 
that the linear equations can in this case be reduced to a well known recurrence 
formula of the Bernoulli numbers. Accordingly, the terms of the corresponding 
Euler-Maclaurin formula are exactly reproduced. 

In Sectional the Gaussian hypergeometric series 2Fi{a,b;c; z) is treated. 
Since the terms of this series depend on three parameters and one argument, a 



4 Ernst Joachim Weniger 



closed form solution of the linear equations seems to be out of reach. Instead, 
approximations are computed symbolically with the help of the computer 
algebra system Maple. The practical usefulness of these approximations is 
demonstrated by some numerical examples. 

In Section [7| the divergent asymptotic inverse power series for the expo- 
nential integral Ei{z) is treated. Again, the linear equations are solved sym- 
bolically with the help of Maple, and the practical usefulness of these solutions 
is demonstrated by some numerical examples. 



2 The Euler-Maclaurin Formula 



The derivation of the Euler-Maclaurin formula is based on the assumption 
that g{x) is a smooth and slowly varying function. Then, J^^ g(x)dx with 
M, iV e Z can be approximated by the finite sum ^g{M) + g{M -I- 1) -I- ■ • • -I- 
g{N — 1) -I- ig(iV). This finite sum can also be interpreted as a trapezoidal 
quadrature rule. In the years between 1730 and 1740, Euler and Maclaurin de- 
rived independently correction terms to this quadrature rule, which ultimately 
yielded what we now call the Euler-Maclaurin formula (see for example |19[ 
Eq. (1.20)]): 



v=M 



gix)dx + -[g{M) + giN)] 



73 



g(2j-i)(iV)- 5(2^-1) (A/) 



Rk{g) 



{2k)\ 



N 



M 



B2k{x^\x\)g'^'"'\x)Ax. 



Rkig), (6a) 
(6b) 



Here, g'^"^\x) is the rn-th derivative, [x\ is the integral part of x, Bm{x) is 
a Bernoulli polynomial defined by the generating function te^*/(6* ~ 1) = 
SJ^o Bn{x)t" /n\, and Bm = -Bm(O) is a Bernoulh number. 

It is not a priori clear whether the integral Rk{g) in ® vanishes as fc ^ 00 
for a given function g{x). Thus, the Euler-Maclaurin formula may lead to an 
asymptotic expansion that ultimately diverges. In this article, it is always 
assumed that the Euler-Maclaurin formula and related expansions are only 
asymptotic in the sense of Poincare. 

Although originally used to express the in the early 18th century still unfa- 
miliar integral in terms more elementary quantities, the Euler-Maclaurin for- 
mula is now often used to approximate the truncation error r„ = — 'Yl^=n+i 
of a slowly convergent monotone series by an integral plus correction terms. 
The power and the usefulness of this approach can be demonstrated convinc- 
ingly via the Dirichlet series for the Riemann zeta function. 
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The terms {ly + l)^** of the Dirichlet series are obviously smooth and 
slowly varying functions of the index i/, and they can be differentiated and 
integrated easily. Thus, the application of the Euler-Maclaurin formula © 
with M = n + 1 and = cx) to the truncation error of the Dirichlet series 
yields: 

k 



Here, (s)m = s{s + 1) ■ ■ ■ {s + m - 1) = r{s + Tn)/r{s) with s 6 C and m e No 
is a Pochhammer symbol. 

In Tables 8.7 and 8.8, p. 380] and in Section 2] it was shown that a 
few terms of the sum in (|7a|l suffice for a convenient and reliable computation 
of with s = 1.1 and s = 1.01, respectively. For these arguments, the 
Dirichlet series for ({s) converges so slowly that it is practically impossible to 
evaluate it by adding up its terms. 

In order to understand better its nature, the Euler-Maclaurin formula © is 
rewritten in a more suggestive form. Let us set M = n + 1 and N = oo, and let 
us also assume liuiN^ao g{N) = limjv^oo = ^^^n^oc g"{N) = • • • = 0. 

With the help of Bq ^ 1. Bi = -1/2, and B2n+i = with n e N (sec for 
example ^1 p. 3]), we obtain: 



di'^) = -Bo / g{x) dx + 

v=n+l "'"+1 n=l 



111 



^m(5) = V4r / B,4x-[x\)g^"'\x)dx, meN. (8b) 

In the same way, we obtain for the Euler-Maclaurin approximation {T)) to the 
truncation error of the Dirichlet series: 

- f: (. + ir = E ^''^'"'^f^-^''^ (n + 2)^— ^ + fi^(n,.), 

i/=n-|-l fi=0 ^' 

(9a) 

(-l)"'(g)„, p B^{x-[x\) ^ 
iim\n,s) = — -7 / — — — -- — -ax, meN. (9b) 

The reformulated Euler-Maclaurin approximation ® looks suspiciously 
like a truncated expansion of the truncation error in terms of the asymptotic 
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sequence {{n + 2)^^}^q of inverse powers. An analogous interpretation of 
the reformulated Euler-Maclaurin formula ^ is possible if we assume that 
the quantities J^_^_i g{x)dx, g{n)^ g'in), g" {n), . . . form an asymptotic sequence 
{0A'('^)}^o as n — > oo according to 

/>oo 

Go{n) = / g{x)dx, (10a) 

Jn+l 

G^.{n) - g'^^-'\n), fieE. (10b) 

The expansion of the truncation error — X^S^n+i di'^) terms of the 
asymptotic sequence {^m("')}^o a-ccording to ^ has the undeniable ad- 
vantage that the expansion coefficients do not depend on the terms g(i^) and 
are explicitly known. The only remaining computational problem is the de- 
termination of the leading elements of the asymptotic sequence {^^(n)}^_|^. 
In the case of the Dirichlet series this is trivially simple. Unfortunately, 
the terms of most series expansions for special functions are (much) more 
complicated than the terms of the Dirichlet series In those less fortunate 
cases, it can be extremely difhcult to do the necessary differentiations and 
integrations. Thus, the construction of the asymptotic sequence {0p("-)}^_o 
may turn out to be an unsurmountable problem. 



3 Asymptotic Approximations To Truncation Errors 

Let us assume that we want to construct an asymptotic expansion of a special 
function f{z) as z — > oo. First, we have to find a suitable asymptotic sequence 
{cpj {z)}'j!Lq. Obviously, {(pj {z)}J?^q must be able to model the essential features 
of f{z) as z ^ oo. On the other hand, {'Pj{z)}'j!^Q should also be sufficiently 
simple in order to facilitate the necessary analytical manipulations. In that 
respect, the most convenient asymptotic sequence is the sequence {z~^}°^q 
of inverse powers, and it is also the one which is used almost exclusively in 
special function theory. An obvious example is the asymptotic series Q. 

The behavior of most special functions as z — > oo is incompatible with an 
expansion in terms of inverse powers. Therefore, an indirect approach has to 
be pursued: For a given /(z), one has to find some g(z) such that the ratio 
f{z)/g{z) admits an asymptotic expansion in terms of inverse powers: 

oo 

f{z)/g{z) ^ ' (11) 

Although /(z) cannot be expanded in terms of inverse powers {z~^°^q, it can 
be expanded in terms of the asymptotic sequence {ff(z)/z^}°^Q. The asymp- 
totic series (jSJ is of the form of lfTT|) with /(z) = Ei{z) and g{z) — exp(— z)/z. 
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It is the central hypothesis ol this article that such an indirect approach 
is useful for the construction of asymptotic approximations to remainders 
r„ = — jy^^n+i of infinite series as n — > oo. Thus, instead of trying to use 
the technically difficult Euler-Maclaurin formula we should try to find 
some pn such that the ratio Tn/ Pn admits an asymptotic expansion as n — > cxd 
in terms of inverse powers {{n + ol)~^}^q with a > 0. 

A natural candidate for /?„ is the first term a„+i neglected in the partial 
sum s„ = X]"=o '^1'' some cases it is better to choose instead pn = an or 

Pn — {n+ajan+i with a > 0. Moreover, the terms a„+i and the remainders r„ 
of an infinite series are connected by the inhomogeneous difference equation 

= - r„ = a„+i , n e Nq . (12) 

In Jagerman's book [HI Chapter 3 and 4], solutions to difference equations of 
that kind are called Norlund sums. 

If we knew how to solve (I12|l efficiently and reliably for essentially arbitrary 
inhomogeneities a„_|_i, all problems related to the evaluation of infinite series 
would in principle be solved. Unfortunately, this is not the case. Nevertheless, 
we can use (|12|l to construct the leading terms of an asymptotic expansion of 
r„/a„+i or of related expressions in terms of inverse powers. 

For that purpose, we make the following ansatz: 

ri"-) = -a„+i V , . , neNo, meN, a>0. (13) 

This ansatz, which is inspired by the theory of converging factors [Tlll3) and by 
a truncation error estimate for Levin's sequence transformation |10| proposed 
by Smith and Ford [TBI Eq. (2.5)] (see also [22 Section 7.3] or [Hj Section 
IV]), is not completely general and has to be modified slightly both in the case 
of the Dirichlet series for the Riemann zeta function, which is discussed 
in Section [3] and in the case of the divergent asymptotic series (0 for the 
exponential integral, which is discussed in Section {7\ Moreover, the ansatz 
lfO|l does not cover the series expansions of all special functions of interest. 
For example, in |23| a power series expansion for the digamma function ip{z) 
was analyzed whose truncation errors cannot be approximated by a truncated 
power series of the type of ((T^ . Nevertheless, the examples considered in this 
article should suffice to convince even a sceptical reader that the ansatz H13f) 
is indeed computationally useful. 

We cannot expect that the ansatz satisfies the the inhomogeneous 
difference equation (fl^ exactly. However, we can choose the unspecified co- 
efficients 7^™^ in H13|l in such a way that only a higher order error remains: 



(m) (m) m m m 

' n+1 _ \ " l/j, (ln+2 \ ^ 



(14) 



l + 0(n-™-'), n^oo. (15) 
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The approach of this article depends crucially on the assumption that the ratio 
O'n+2/ dn+i can be expressed as an (asymptotic) power series in l/(n + a). If 
this is the case, then the right-hand side of H14|) can be expanded in powers 
of l/(n + a) and we obtain: 

„(™) _ ™ r^-m) 

= > 7 r hOn , n —f 00 . (16) 

a„+i ('^ + 

Now, (|15|l implies that we have solve the following system of linear equations: 

C(") =5^0, 0</i<m. (17) 

Since C^™^ with < /i < m contains only the unspecified coefficients 7q'"\ . . . , 
7}™'' but not 7^™i • ■ ■ , 7m"'' , the linear system (|17ll has a triangular structure 
and the unspecified coefficients 70™'' , • ■ ■ , 7^™'' can be determined by solving 
successively the equations Cq™' = 1, cj™'' =0, . . . , cln'' = 0. 

Another important aspect is that the linear equations (|17|) do not depend 
explicitly on m, which implies that the coefficients 7^™'' in H13|) also do not 
depend explicitly on m. Accordingly, the superscript m of both C^™'' and 7/,™^ 
is superfluous and will be dropped in the following Sections. 



4 Numerical Analytic Continuation 

Divergent asymptotic series of the type of (|ll|l can be extremely useful com- 
putationally: For sufficiently large arguments z, truncated expansions of that 
kind are able to provide (very) accurate approximations to the corresponding 
special functions, in particular if the series is truncated in the vicinity of the 
minimal term. If, however, the argument z is small, truncated expansions of 
that kind produce only relatively poor or even completely nonsensical results. 

We can expect that our asymptotic expansions in powers of l/(n-t-Q;) have 
similar properties. Thus, we can be confident that they produce (very) good re- 
sults for sufficiently large indices n, but it would be overly optimistic to assume 
that these expressions necessarily produce good results in the nonasymptotic 
regime of moderately large or even small indices n. 

Asymptotic approximants can often be constructed (much) more easily 
than other approximants that are valid in a wider domain. Thus, it is desirable 
to use asymptotic approximants also outside the asymptotic domain. This 
means that we would like to use our asymptotic approximants also for small 
indices n in order to avoid the computationally problematic asymptotic regime 
of large indices. Obviously, this is intrinsically contradictory. We also must 
find a way of extracting additional information from the terms of a truncated 
divergent inverse power series expansion. 
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Often, this can be accomplished at low computational costa by converting 
an inverse power series X^^o*^"/-^" ^ factorial series Y^^=o'^n/ {z)n- Fac- 
torial series, which had already been known to Stirling POl p. 6], frequently 
have superior convergence properties. An example is the incomplete gamma 
function r(a, z), which possesses a divergent asymptotic series of the type of 
l|TT]l Eq. (6) on p. 135] and also a convergent factorial series [3 Eq. (1) 
on p. 139]. Accordingly, the otherwise so convenient inverse powers are not 
necessarily the computationally most effective asymptotic sequence. 

The transformation of an inverse power series to a factorial series can be 
accomplished with the help of the Stirling numbers of the first kind which 
are normally defined via the expansion (z — n + 1)„ ~ X)"=o ^'■^^("■j ^)z^ of 
a Pochhammer symbol in terms of powers. As already known to Stirling (see 
for example 20, p. 29] or TM. Eq. (6) on p. 78]), the Stirling numbers of the 
first kind occur also in the factorial series expansion of an inverse power: 

1 ^ (-l)«S(l)(fc + K,/c) , 

This infinite generating function can also be derived by exploiting the well 
known recurrence relationships of the Stirling numbers. 

With the help of (|18|) , the following transformation formula can be derived 
easily: 



r!=0 ^ k=2 



(^)l ^2 (^)fc t^l 



Let us now assume that the coefficients 7^ with < < to of a truncated 
expansion of r„/a„_|_i in powers of 1/(?t. + a) according to (|13|1 are known. 
Then, H19I) implies that wc can use the transformation scheme 




^ = 0,1, 

^(-1)^+^sW(a.-1,^^-1)7., /i>2, 
. 1^=1 

to obtain instead of (|13|l the truncated factorial series 



(20) 



f*^™' = -a _Li 

' n "n+l 



.fj.=0 



n^oo. (21) 



Pade approximants, which convert the partial sums of a formal power series 
to a doubly indexed sequence of rational functions, can also be quite helpful. 
They are now used almost routinely in applied mathematics and theoretical 
physics to overcome convergence problems with power series (see for example 
the monograph by Baker and Graves-Morris P] and references therein). The 
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ansatz (|13|) produces a truncated series expansion of r„/a„_|_i in powers of 
l/(n + a), which can be converted to a Pade approximant, i.e., to a rational 
function in l/{n + a). 

The numerical results presented in Sections and [7| that the conversion to 
factorial series and Pade approximants improves the accuracy of our approx- 
imants, in particular for small indices n. 



5 The Dirichlet Series for the Riemann Zeta Function 

In this Section, an asymptotic approximation to the truncation error of the 
Dirichlet series for the Riemann zeta function is constructed by suitably adapt- 
ing the approach described in Sectional In the case of the Dirichlet series 
we have: 

n 

s„ = ^(,. + 1)-% (22) 

oc oc / \ —s 

r„ = - ^ (i. + l)- = -(n + 2)-^ 1 + -^) , (23) 

zir„ = (^^ + 2)-^ (24) 

It is an obvious idea to express the infinite series on the right-hand side 
of H23|l as a power series in l/{n -I- 2). If z/ < n -I- 2, we can use the binomial 
series (1 -t- 2)" = iFo(-a; -z) = J2^^o (^)^"' 112 P- 38], which converges for 
\z\ < 1. We thus obtain 



[1 + Win + 2)]- ^ £ IfM [_^/(,, + 2)r . (25) 



{s 

rn—0 



The infinite series converges if j//(n + 2) < 1. Thus, an expansion of the right- 
hand side of (|23|l in powers of l/(n -I- 2) can only be asymptotic as n — > 00. 
Nevertheless, a suitably truncated expansion suffices for our purposes. 

In the case of the Dirichlet series for the Riemann zeta function, we cannot 
use ansatz H13II . This follows at once from the relationship 

Z\n" = (n+l)"-n" = an°'-^ + 0{n°'^^) , n -> oo . (26) 

Thus, we make the following ansatz, which takes into account the specific 
features of the Dirichlet series : 



ri") - (n + 2)i-- V men, neNo- (27) 

■^-^ (n + 2)** 



This ansatz is inspired by the truncation error estimate for Levin's u trans- 
formation (see also EH Section 7.3] or "I^, Section IV]). 
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As in Section 121 the unspecified coefficients 7^ are cliosen in sucli a way 
tfiat oniy a higlier order error remains: 



For that purpose, we write: 



n 00 . 



'n+l 

(n + 2)i"« 



n + 3 



71 + 2 



1 — s m 



E 



7m 



E 



7m 



{1 - [l + l/(n + 2)]i-^-^} 



Witii the heip of the binomiai series |12[ p. 38], we obtain: 

(s + A*-1)a+i (-1)^ 



1- [l + l/(n + 2)] 



1 — S — /J, 



E 

A=0 



(A + 1)! (n + 2)^+1 



Inserting into (|Sn|l yieids: 



E 

M=0 



+ + 2)] — } = 



(s + M-1)a+i (-1)' 



(n + 2)A' ^ (A + 1)! (n + 2)-^+i 



min(z^,rn) 



E(-+2)— ^ E 



(1 - 5 - i^)A+l7i^-A 

(A + 1)! 



Thus, we obtain the following truncated asymptotic expansion: 



„(™) _ 
(n + 2)-''- 



E 

A=0 



-^(n + 2)-'' 

M=0 

(1 - S - Ai)A+l7M-A 



(A + 1)! 



0( 



(28) 

(29) 
(30) 

(31) 



(32) 
(33) 



(34) 



The unspecified coefficients 7^ have to be determined by solving the following 
system of linear equations, whose triangular structure is obvious: 



E 



(1 -s-/z)a+i7m-a 



< /Lt < m . 



(35) 



A=0 



(A + 1)! 

For a more detailed analysis of the linear system 1)351) . let us define via 

M e No . (36) 



7. = (-1)''^/?., 
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Inserting (|3(j|l into (|35|l yields: 

(A + 1)! O^^A)^ = °' ^ ^ 

A — 

Next, we use (1 — s — /i)A+i = (— l)^^^(s)p and obtain 

A. (A + 1)! - ( 1^ (^)mI. (38) 

- (^ + 1)! A. (^_a + l)w'^'^ ^^^^ 

Thus, the hnear system ()35|l is equivalent to the well known recurrence formula 

= 0, neN, (41) 



of the Bernoulli numbers (see for example |19[ Eq. (1.11)]) together with the 
initial condition Bq — 1. Thus, the ansatz (|27|l reproduces the finite sum 17a|l 
or (|9a|l of the Eulcr-Maclaurin formula for the truncation error of the Dirichlet 
series, which is not really surprising since asymptotic series are unique, if they 
exist. Only the integral ITbp or IjQbll ) cannot be reproduced in this way. 



6 The Gaussian Hypergeometric Series 

The simplicity of the terms of the Dirichlet series ^ facilitates the deriva- 
tion of explicit asymptotic approximations to truncation errors by solving a 
system of linear equations in closed form. A much more demanding test for 
the feasibility of the new formalism is the Gaussian hypergeometric series ^ , 
which depends on three parameters a, &, and c, and one argument z. 

Due to the complexity of the terms of the Gaussian hypergeometric series 
(O, there is little hope in obtaining explicit analytical solutions to the linear 
equations. From a pragmatist's point of view, it is therefore recommcndable 
to use computer algebra systems like Maple and Mathematica and let the 
computer do the work. 

In the case of a nonterminating Gaussian hypergeometric series, we have: 
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i^— n+1 ^ ' 



(a)„+i(b)„+i ^ (a + n + l)^(6 + n + l)„ ^ 



(c)„+i(n + l)! {c + n + l)^{n + 2)^ 

{a)n+l{b)r 



(c)„+i(n + 1)! 

Since [{a + n+ l)„{b + n+ l)^]/[{c + n+ l)„{n + 2)^] can be expressed as 
a power series in l/(n + 1), the following ansatz make sense: 



r, 



(m 



(a)„+i(6)n+i +1 -A 7m 



(c)„+i(n+l)! 

TO e N, n e No , |z| < 1 . (45) 

Again, we choose the unspecified coefficients 7^ in (|45|l in such a way that 
only a higher order error error remains: 



(a)„+i(6)„ 



ii- [1 + 0(n-"-i)] , n ^ 00 . (46) 



(c)„+i(n + l)! 

This convergence condition can be reformulated as follows: 



[(a)„+i(&)„+i^"+i]/[(c)„+i(n+l)!] 

7m _ (a + n + l)(b + ?t + l) ^ >A 7^. /^^x 
Z^(n + 1)^ (c + n+l)(n + 2) ^ (n + 2)^ ^ 

= 1 + 0(n-"-i), n^c5o. (48) 

Now, we only have to do an asymptotic expansion of (|47l) in terms of the 
asymptotic sequence {l/(n + 1)-'}°^q as n ^ 00. This yields: 



[(a)„+i(6)„+iz"+i]/[(c)„+i(n + l)!] 

- E7;^ + 0(n— 1), n-.oo. (49) 

We then obtain the following system of coupled linear equations in the un- 
specified coefficient 7^ with < /i < to: 
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^ S^o, 0<fi<m. (50) 

As discussed in Section |3| a coefficient with < /i < m contains 
only the unspecified coefiicients 70, . . . , 7^ but not 7^+1 . . . , 7m- Thus, the 
symbofic solution of these linear equations for a Gaussian hypergeometric 
function 2-^1 (a, b; c; z) with unspecified parameters a, b, and c and unspecified 
argument z is not particularly difficult for a computer algebra system, since 
the unspecified coefficient 7^ can be determined successively. The following 
linear equations were constructed with the help of Maple 8: 

Co - (1-2)70 = 1, (51a) 
Ci = (c-a-&+l)z7o + (l-z)7i = 0, (51b) 

C2 = [ic-b+l)a+{c+l)b-l-c-c^]zjo 

+ (c + 2-6-a)z7i + (l-z)72 = 0, (51c) 

C3 = {[{c+l)b-l-c^c'^]a-{l + c + c^)b + c^ + c^ + c+l}zjo 
+ [(c + 2-6)a+(c + 2)6-3-c2-2c]z7i 
+ (3-5-a + c)z72 + (l-z)73 = 0. (51d) 

This example shows that the complexity of the coefiicients in (|49|l increases 
so rapidly with increasing index fj, that a solution of the linear equations 
H5U|I becomes soon unmanageable for humans. This is also confirmed by the 
following solutions of H51|l obtained symbolically with the help of Maple 8: 



(2-1)70 = 1, (52a) 

(z_l)2^-^ = z(a + 6-c-l), (52b) 
(2-1)^72 = z{[a^ + (b-c-2)a + b^ - {c + 2)b+l + 2c]z 

+ {b-c-l)a-{c+l)b+l + c + c^}, (52c) 



(z - 1)S3 = z {[a^ + (6 - c - 3) + (6^ - (c + 3) 6 + 3 + 3 c) a 
+ 6''' - (c + 3) 6^ + (3 + 3 c) 6 - 1 - 3 c] 
+ [(2 - 2 c - 3) + (2 6^ _ (4 c + 8) 6 + 2 + 7 + 8 c) a 
- (2 c + 3) 6^ + (2 + 7 + 8 c) 6 - 4 - 5 - 7 c] z 
+ [-{c+l)b+l + c^ + c] a + {l + c^ + c)b 
-\-(?~c-e'}. (52d) 

The solutions (|52|l . which are rational in z, demonstrate quite clearly a 
principal weakness of symbolic computing. Typically, the results are compli- 
cated and poorly structured algebraic expressions, and it is normally very 
difficult to gain further insight from them. Nevertheless, symbolic solutions of 
the linear equations H5()|l are computationally very useful. 

For the Gaussian hypergeometric series with a = 1/3, b — 7/5, c — 9/2, 
and z — —0.85, Maple 8 produced for m = 8 and n ~ 1 the following results: 
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ri = -0.016 412 471, (53a) 

02 [4/4] = -0.016 410 482, (53b) 

02 ff' = -0.016 414 203, (53c) 

02 ^ = -0.004 008 195 . (53d) 



It is in my opinion quite remarkable that for n — 1, which is very far away 
from the asymptotic regime, at least the Pade approximant 02 [4/4] and the 

truncated factorial series a2f\ agree remarkably well with the "exact" trun- 

(s) 

cation error ri. In contrast, the truncated inverse power series a2r\ produces 
a relatively poor result. For n — 10, which possibly already belongs to the 
asymptotic regime. Maple 8 produced the following results: 

rio = 0.000 031 925 482 , (54a) 
an [4/4] = 0.000 031 925 482 , (54b) 

an f = 0.000 031 925 483 , (54c) 

an rJo^ = 0.000 031 925 471 . (54d) 

Finally, let me emphasize that the formalism of this article is not limited 
to a Gaussian hypergeometric series (0), but works just as well in the case of 
a generalized hypergeometric series p+iFp (ai , . . . , ap+i , . . . , /3p; z) . 



7 The Asymptotic Series for the Exponential Integral 

The divergent asymptotic series (EJ for the exponential integral -Ei(z) is prob- 
ably the most simple model for many other factorially divergent asymptotic 
inverse power series occurring in special function theory. Well known examples 
are the asymptotic series for the modified Bessel function K^{z), the comple- 
mentary error function erfc(z), the incomplete gamma function r{a,z), or 
the Whittaker function W,^^p{z). Moreover, factorial divergence is also the 
rule rather than the exception among the perturbation expansions of quan- 
tum physics (see for a condensed review of the relevant literature) . 

The exponential integral Ei{z) can also be expressed as a Stieltjes integral: 

r°° e"*di 

If z < 0, this integral has to be interpreted as a principal value integral. 

In the case of a factorially divergent inverse power series, it is of little 
use to represent the truncation error r„(z) by a power series as in H43II . If, 
however, we use X]"=o s;"^ = [1 — a;"+^]/[l — x] in H55|l . we immediately obtain: 
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n 

Sn{z) = Y,{-l/zriy\, (56) 

Arniz) = (-l/zr+i(n+l)!. (58) 

Because of the factorial growth of the coefficients in (|3Jl, it is advantageous to 
use instead of H13|l the following ansatz: 

m 

r(")(z) = (-l/z)"n! ^ mSN, n G No . (59) 

Again, we choose the unspecified coefficients 7^ in (|59|l in such a way that 
only a higher order error error remains: 

Z\ri™)(z) = (-f/z)"+i(n + l)! [1 + 0(n-™-i)] , n -> 00 . (60) 

This convergence condition can be reformulated as follows: 

(m) / \ (m) / \ m m 

r'„+[{z)-r'n '{z) ^ -z 7p _ 7m /g.N 

(-l/z)"+i(n + l)! n + \^{n^\Y ^{n^2Y ^ ' 

= 1 + 0{n-"'-^) , n ^ 00. (62) 

Next, we do an asymptotic expansion of the right-hand side of (|61|l in terms 
of the asymptotic sequence {l/(ri + l)-'}|^o ^ ~^ "^^^^ yields: 

— — -r, r: = > -, h 0(n I, n ^ 00 . (63) 

(-l/z)"+i(n-t- 1)! ^("+1)^ 

Again, we have to solve the following system of linear equations: 

= <5^o , 0<^i<m. (64) 

The following linear equations were constructed with the help of Maple 8: 

Co = -70 = 1 , (65a) 

Ci = -Z70-71 = 0, (65b) 

C2 - (1 - ^) 71 - 72 = , (65c) 

C3 = -71 + (2 - ;2) 72 - 73 = 0, (65d) 

C4 = 71 - 372 + (3 - z) 73 - 74 = . (65e) 

If we compare the complexity of the equations (|65|l with those of H51|l . we 
see that in the case of the asymptotic series © for the exponential integral 
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there may be a chance of finding exphcit expressions for the coefficients 7^. 
At least, the solutions of the linear system Ht)5|l obtained symbolically with 
the help of Maple 8 look comparatively simple: 

70 = -1, (66a) 

-fi — z , (66b) 

72 = -{z - 1) z , (66c) 

73 = (z2-3z + l)z, (66d) 

74 = ~{z^ -6z^ + 7z-l)z. (66e) 

Of course, this requires further investigations. 

The relative simplicity of the coefficients in (|66|l offers other perspectives. 
For example, the [2/2] Pade approximant to the truncated power series in 
H59|l is compact enough to be printed without problems: 

3 - z 2 
-1 + 



^ -+1 jn + l)^ ^ n -n + zn + z 

^'^ 2z-3 z^-2z + 2 n^-n + 2zn + z^' ^ ' 

n+1 (ri + l)2 

Pade approximants to the truncated power series in H59() seem to be a new 
class of approximants that are rational in both n and z. 

For the asymptotic scries (jSJ for E\{z) with z = 5, Maple 8 produced for 
m = 16 and n = 2 the following results: 

r2 = 0.027 889, (68a) 
02 [8/8] = 0.027 965, (68b) 

02 4^^^ = 0.028 358 , (68c) 

02^^^^ = -177.788. (68d) 

The Pade approximant 02 [8/8] and the truncated factorial series a^r'^^'' agree 
well with the "exact" truncation error r2, but the truncated inverse power 
series 02^2^^'' is way off. For n = 10, all results agree reasonably well 

rio = 0.250 470 879 , (69a) 
aio[8/8] = 0.250 470 882, (69b) 

aiofJo®^ = 0.250 470 902, (69c) 

aiorJo^^ = 0.250 470 221. (69d) 



8 Conclusions and Outlook 

A new formalism is proposed that permits the construction of asymptotic 
approximations to truncation errors — — X^l^n+i of infinite series for 
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special functions by solving a system of linear equations. Approximations to 
truncation errors of monotone series can be obtained via the Euler-Maclaurin 
formula. The formalism proposed here is, however, based on different assump- 
tions and can be applied even if the terms of the series have a comparatively 
complicated structure. In addition, the new formalism works also in the case 
of alternating and even divergent series. 

Structurally, the asymptotic approximations of this article resemble the 
asymptotic inverse power series for special functions as z ^ oo, since they are 
not expansions of r„, but rather expansions of ratios like rn/a,n+i, fn/o-m or 
r-a /[{P'+a)an+i] with a > 0. This is consequential, because it makes it possible 
to use the convenient asymptotic sequence {l/(n + q)-'}°^q of inverse powers. 
This greatly facilitate the necessary analytical manipulations and ultimately 
leads to comparatively simple systems of linear equations. 

As shown in Section |S1 the new formalism reproduces in the case of the 
Dirichlet series for the Riemann zeta function the expressions (|7all or (|9a|l 
that follow from the Euler-Maclaurin formula. The linear equations (|35|l are 
equivalent to the recurrence formula 1)41(1 of the Bernoulli numbers. Thus, only 
the integral ((7bp or ((9bp ) cannot be obtained in this way. 

Much more demanding is the Gaussian hypergeometric series l|2Il, which 
is discussed in Sectional The terms of this series depend on three in general 
complex parameters a, 6, and c and one argument z. Accordingly, there is little 
hope that we might succeed in finding an explicit solution to the linear equa- 
tions. However, all linear equations considered in this article have a triangular 
structure. Consequently, it is relatively easy to construct solutions symboli- 
cally with the help of a computer algebra system like Maple. The numerical 
results presented in Section also indicate that the formalism proposed in 
this article is indeed computationally useful. 

As a further example, the divergent asymptotic series Q for the exponen- 
tial integral Ei(z) is considered in Section [3 The linear equations are again 
solved symbolically by Maple. Numerical results are also presented. This ex- 
ample is important since it shows that the new formalism works also in the 
case of factorially divergent series. The Euler-Maclaurin formula can only 
handle convergent monotone series. 

Although the preliminary results look encouraging, a definite assessment of 
the usefulness of the new formalism for the computation of special functions 
is not yet possible. This requires much more data. Consequently, the new 
formalism should be be applied to other series expansions for special functions 
and the performance of the resulting approximations should be analyzed and 
compared with other computational approaches. 

I suspect that in most cases it will be necessary to solve the linear equa- 
tions symbolically with the help of a computer algebra system like Maple. 
Nevertheless, it cannot be ruled out that at least for some special functions 
with sufficiently simple series expansions explicit analytical solutions to the 
linear equations can be found. 
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Effective numerical analytic continuation methods are of considerable rel- 
evance for the new formalism which produces asymptotic approximations. We 
cannot tacitly assume that these approximations provide good results outside 
the asymptotic regime, although it would be highly desirable to use them also 
for small indices. In Section 0] only factorial series and Pade approximants 
are mentioned, although many other numerical techniques are known that 
can accomplish such an analytic continuation. Good candidates are sequence 
transformations which are often more effective than the better known Pade 
approximants. Details can be found in books by Brezinski and Redivo Zaglia 

Sidi 22], or Wimp |26|, or in a review by the present author P^. 
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